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ABSTRACT 



The performance of Loran-C receivers is severely limited by 
contamination introduced by multipath propagation effects. Echo removal 
by discrete generalized linear filtering is considered as a method of 
reducing the effects of the contamination. The signals can be described 
as consisting of repeated, known, composite signals with short echo 
arrival times, in noise. Since the signals are repeated, and stable, 
they can be averaged to increase the effective signal to noise ratio. 

The homomorphic deconvolution technique is applied to simulated signals 
and a determination of the averaging tine needed for satisfactory 
performance is made. 

Aside from the Loran-C application, the report presents a detailed 
study of the computational aspects involved in deriving the complex 
cepstrum. Since the signal is known, the effects of noise and aliasing 
can be precisely identified. In particular, conventional phase 
unwrapping methods are questioned for certain classes of signals. 
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Chapter 1 
INTRODUCTION 

The digital signal processing technique referred to a Homomorphic 
Deconvolution has met with some success, in a number a diverse appli- 
cations, in suppressing signal contamination introduced by a multi- 
path propagation medium. Another possible application is to Loran-C, 
a low frequency electronic navigation system, in which the performance 
is severely limited by the effects of such contamination. It is the 
purpose of this thesis to study the applicability of the technique to 
various aspects of Loran-C signal processing. 

In this Chapter, a brief description of the Loran-C system and 
Homomorphic Deconvolution is presented. In Chapter 2, the task of the 
homomorphic filter in the Loran-C case is formulated. In Chapter 3 > 
the filtering technique is developed and the results of simulated 
filtering in the noiseless case are presented. Chapter 4 presents the 
theory and a simulated example of the effects of signal distortion. 
Chapter 5 includes a theory of the effect of white gaussian noise on 
the processing and includes examples to verify the theory. In 
Chapter 6, the effects of the various assumptions and simplifications 
made are discussed and the results of the preceeding chapters are 
examined to determine the applicability of the method. The examples 
presented are simulations performed on an IBM 370 general purpose 



computer 



1.1 Loran-C System Description 
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Loran-C is an electronic navigation system which belongs to a 
broad class of systems referred to as hyperbolic, a characterization 
based on straightforward geometric considerations. If radio signals 
are simultaneously transmitted from two different locations, they will, 
in general, arrive at some arbitrary third location at two different 
times. To the extent that the speeds of propagation along the two 
paths are equal and known, a measurement of the difference in times 
of arrival provides a corresponding measure of the difference in 
distances involved. The locus of points with a constant difference 
in distance from two reference points is a hyperbola. A receiver 
which determines this time difference and identifies the proper 
reference points provides information sufficient to define a navi- 
gational line of position. Tito or more such lines of position can 
be used to obtain a fix - the location of the receiver. 

Various hyperbolic systems with origins in the late 1930's have 
evolved into today's alternatives. Of these, it may be said, the one 
which meets the most general operational requirements is the Loran-C 
system. The advantages of Loran-C can, in short, be attributed to 
the following. 

1. Use of a frequency band (90 to 110 kHz) which provides 
for extended range and exhibits excellent propagation 
properties. 

. Use of a pulsed signal format. This increases average 
transmitted power and permits signal coding to facilitate 



2 
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transmitting station identification and to reduce the effects 
of interference, 

3. Improved measurement precision obtained by coarse measure- 
ment of the envelope followed by fine tune measurements based 
on the signal rf carrier.^ 

In Loran-C, the hyperbolic navigation concept is realized in a 
somewhat complex manner. There are a number of modifications to the 
simple scheme as previously described. While the modifications make the 
system appear involved, they are seen to be well founded and straight- 
forward when systematically examined. 

A first consideration is the operation of the transmitting stations. 
These are grouped into chains whose configurations, subject to practical 
constraints, provide geometrically optimal accuracy and coverage. An 
example of a simple chain to analyze is the so-called triad which 
consists of stations denoted master, slave X and slave Y. The trans- 
missions consist of groups of pulses on a 100 kHz carrier frequency: 
the master transmits nine pulses in a group while the slaves each 
transmit eight pulses in a group. In practice, the stations do not 
transmit simultaneously. The slave stations incorporate known trans- 
mission delays relative to the master so that, within the service 
area of a triad, reception is typically as depicted in figure 1-1. 
Transmission timing is controlled so that the groups are received in the 
order and intervals shown. The transmission delays guarantee that there 
will be no overlap of the groups. In this example, the transmission 
sequence is repeated every 30,000 microseconds - the chain pulse 
group repitition interval. Different chains are distinguished by 
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50,000 psec 



Chain Pulse Group Interval 



Master 
Pulse Group 
Interval 



Slave X 
Pulse Group 
Interval 



Slave Y 
Pulse Group 
Interval 



-A h- 

1000 |isec 



Figure 1-1 Loran-C Chain Fulse Transmission Scheme 



different group repitition intervals. Consequently, by means of a 
correlation type operation applied over a number of repitition intervals, 
a receiver can discriminate between signals transmitted from adjacent 
chains. 

The stations also employ a pulse phase code in the transmission. 

This allows enhanced discrimination between signals from competing 
chains and overcomes some of the problems introduced by skywave 
contamination. 

The term skywave refers to the multipath propagation phenomenon 
depicted in figure 1-2. The high accuracy and coverage advantages of 
Loran-C depend on the use of measurements obtained from signals 
arriving via the stable groundwave propagation mode at extended ranges. 
The reception of skywaves, whose arrival times are highly dependent 
upon unpredictable ionospheric properties, causes severe signal 
processing problems. The nature of the important properties of the 
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multipath medium, at Loran-C frequencies, is illustrated in figures 
1-3 and 1-4, It can be seen that the times of arrival can be as short 
as about 35 microseconds after the arrival of the groundwave for strong 
first hop skywaves at extended ranges. Additionally, arrival times can 
be delayed as long as about 1500 microseconds for higher order hops 
under some conditions. An explanation of how these effects are 
overcome requires a further description of the transmitted signal 
format and a brief description of receiver detection methods. 

All transmitting stations attempt to generate pulses of the 

2 -at 

same form which can be described as t e with a such that a peak 
occurs 65 microseconds after the beginning of the pulse 0 Such a pulse 
has an effective duration of 250 to 300 microseconds. An example 
which combines the two most deleterious effects of the multipath 
contamination is illustrated in figure 1-5. For clarity, the magnitudes 




Figure 1-3 Skywave Arrival Times 




Figure 1-4 Skywave Magnitudes 
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of the separate components of the signal have been superimposed. 

The effects of the contamination can be grouped into two 
categories. The groundwaves are overlapped by skywaves associated with 
the same transmitted pulse. These skywaves correspond to low order 
hops and may be many times larger than the groundwave. The groundwaves 
are also overlapped by much delayed, high order hops associated with 
the preceeding transmitted pulse. These skywaves will be fairly 
well attenuated and smaller than the groundwaves they overlap. 

The problems of long delayed skywave interference is easily 
overcome by use of a pulse phase coding. The phase of the trans- 
mitted signal carrier is systematically varied in 180° steps from 
pulse to pulse in accordance with the scheme shown in figure 1-6, 
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Repitition 

Interval 


Any Master Station 
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+ + +'+ + + 


+ - + -+ + 


Second 


+ + + - + - 


+ + + + + + 


Third 


same as first 


same as first 


Fourth 


same as second - 


same as second 


etc* 


etc, 


etc.. 



Pi gure 1-6 Pulse Phase Code 



As an example of the use of the code, consider the effect of the 
sixth hop skywave in figure 1-5 which will cause problems by over- 
lapping the next groundwave. It can be verified from figure 1-6 that 
the code can be represented as a signal, s(t), with the property that 
s(t) and s(t-T) are orthogonal when averaged over two pulse group 
repitition intervals, where T is a multiple of 1000 microseconds. 
Since the effect of the sixth hop skywave is to create a signal 
s(t) +A s(t-T), if the received signal is correlated with s(t), this 
type of interference is removed. 

The effects of the interference caused by low order skywaves 
are not so easily overcome. Use is made of the fact that the 
groundwave travels the most direct path and is therefore guaranteed 
to arrive before any of the skywaves associated with the same trans- 
mitted pulse 0 Consequently the first portion of the groundwave is 
free of contamination of this type and can be used to make the time 
measurements. This portion, however, can not be assumed, a priori, 
to be much more than the first 30 microseconds of the pulse. 

The measurement technique proceeds in two steps. First it is 





Figure 1-7 Envelope Deriver and RF Samplers 
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attempted to identify some point on the leading edge of the pulse. 
Since this point must be within the first JO microseconds after the 
beginning of the pulse, it is clear that this is a severe constraint 
on a pulse which has an effective bandwidth of only 10 kHz. 

A typical technique employed is to process the pulse envelope 
so as to produce a zero crossing at the 25 microsecond point. This 
is accomplished by delaying the received signal 5 microseconds and 
scaling it to produce what is represented by the dashed line in 
figure l-7(a). This is then subtracted from the original signal - 
the solid line in figure l-7(a) to produce the derived envelope 
of figure l-7(b). The zero crossing is then detected and sampling 
strobes are centered at the estimated 25 microsecond point of the 
rf signal as shown in figure 1-7 (c). 

With the strobes located as shown, the middle one will obtain 
a zero sample value while the other two will obtain values with a 
known ratio. If the strobes were not centered properly, different 
values would be obtained, producing an error signal which can be 
used to reposition the strobes. 

In practice, of course, the signals are noisy so that both 
steps in the procedure must be accomplished via measurements on a 
number of pulses - as many as 6k or 128 or more - depending on the 
signal to noise ratio. Typically, the most difficult step in the 
procedure is the envelope zero crossing deriving. The zero crossing 
measurement must be accurate to within 5 microseconds or the fine 
tune measurement that follows will be made on the wrong cycle of 
the carrier. In this regard, a final characteristic of the signals 
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should be presented. 

The propagation paths involved do constitute dispersive media 
so that the signals received are not quite identical to those trans- 
mitted. The distortion in the shape of the envelope is slight - so 
much so that even though there are different propagation paths 
involved, the skywaves may be modelled as perfect echoes of the 
groundwave. The major effect of the propagation paths is the intro- 
duction of a phase shift of the pulse envelope relative to the carrier. 

In effect, the envelope and carrier propagate at different velocities. 
Since the final time difference measurement is based on the carrier 
times, the carrier propagation times are calibrated for measuring 
distances. While the phase shift, for the groundwave, is typically 
not enough, by itself, to cause a cycle ambiguity in the detection, 
it does require that the envelope deriver have accuracy somewhat 
better than 5 microseconds. 

With the preceeding as background, some aspects of system per- 
formance can be presented. Inaccuracies can be introduced by any 
number of sources in the system. These can be loosely grouped into 
a few categories and statistically described. Receiver error will 
have a standard deviation of from 20 to 100 nanoseconds depending on 
the degree of sophistication. Chain synchronization includes the 
results of many effects and may be described as introducing an 
error with a standard deviation of about 50 nanoseconds. Noise 
induced errors can have a standard deviation of from about 20 nano- 
seconds for 0 db SNR to about 250 nanoseconds for -20 db SNR. The 
implications of these times can be seen from an example. If a total 
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error of about 100 nanoseconds is present, this corresponds, via a 
speed of propagation conversion, to a distance error of about 30 
meters. 

The numbers presented above for the error statistics were 
purposely given without mention of important descriptive parameters. 

For example, a figure for the error caused by atmospheric noise is 
completely meaningless without a detailed description of the 
particular receiver detection procedure - particularly the sample 
averaging or integration time involved (50 seconds in the case of the 
above figures). The geometry of the chain and the receiver position 
must also be specified: the fix obtained from lines of position inter- 
secting at right angles can be expected to have an error almost an 
order of magnitude smaller than if the lines were nearly asymptotic 
and at extreme ranges. Without going into great detail, the per- 
formance can be simply summarized, 

1. Position errors vary from a few tens of meters to not 
more than about 200 meters within the coverage region. 

2. The only hyperbolic navigation systems that are competi- 
tive in performance have a coverage range limit of less 
than 200 miles. Loran-C has a nominal coverage range 

of about 1000 miles. 

3. The receiver generally contributes the smallest 
percentage of the total error. 

In spite of the third property, improved receiver performance is 
constantly strived for since it indirectly effects many of the 
other sources of error. Chain synchronization is a single example of 
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such a source. Transmission timing is controlled by passive 
Cesium Beam frequency standards located at each station. As precise 
as these may be, there will be, over a period of time, relative 
phase drifts between the stations of a chain. These phase drifts 
are observed by monitoring stations which use the most sophisticated 
receivers and periodic corrections are ordered to maintain syn- 
chronization to within 100 nanoseconds. Clearly, if more sensitive 
receivers were available, this monitoring process could be made 
more responsive to oscillator phase discrepancies. 

Speculation on means of achieving improved processing schemes 
is tempting and many faceted. One aspect to be considered is that 
an optimum receiver need not totally restrict the measurements to be 
made in the portion of the groundwave that is perfectly free from 
skywave contamination. If, for example, the sampling strobes are 
centered at the 30 microsecond point rather than at the 25 micro- 
second point, the signal to atmospheric noise ratio would be im- 
proved, In such a case, it might result that the third strobe 
would be sampling a signal contaminated by a skywave. In many- cases 
however, the skywave would be of such a low power level at the 
sampling point that the increase in total noise power - atmospheric 
noise plus skywave interference - is more than offset by the in- 
crease in signal power. The implementation of the resulting 
adaptive time filter would require the knowledge of skywave relative 
magnitude and arrival time so that a decision could be made whether 
or not to re-position the sampling points, 

A less optimal, but still improved, receiver would make 
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measurements on the skywave free portion of the signal - but at 
later times in the many cases in which skywave arrival times are 
much more than 30 microseconds after the beginning of the ground wave. 

In this case, the receiver requires only the knowledge of the 
skywave arrival time. More ambitious attempts have been made to 
implement a total power receiver. In this receiver, an attempt is made 
to digitally implement a scheme to measure skywave arrival times, 
magnitudes and phases, and then to subtract them from the received 
signal thereby greatly increasing the signal to noise ratio. 

Another improvement that can be thought of would be allowed by 
a receiver capable of making precise measurements on the envelope of 
the groundwave. This would overcome the problems introduced by the 
different propagation speeds which might cause cycle ambiguity in the 
measurement. Even if a measurement based on the envelope might 
not have quite the precision of the rf measurement and therefore 
not be of use in a navigation receiver, it would be useful in 
calibrating regions in which the envelope - carrier phase shift 
is significant. 

The investigation of an alternate method of accomplishing some 
of these goals is the subject of this report. The method - homo- 
morphic deconvolution - will now be presented. 

1.2 Homomorphic Deconvolution 

In 19 65, A.V. Oppenheim presented a theory of generalized super- 
position which could be applied to certain classes of non-linear 
2 

systems. A simplified description of the theory can be given via 
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an example. Consider processing a signal consisting of two or more 
components which have been combined in some non-linear fashion. By 
means of certain elementary operations, it may be possible to transform 
the non-linear combining operation into a linear one. If the components, 
when so transformed, show a certain distinctiveness, it may be possible 
to process them by conventional linear techniques. One of the components 
can be enhanced relative to the others, or the components can be sepa- 
rated. If the inverse transformation can be realized, the recovered 
signal has been obtained by what is referred to as generalized linear 
filtering. 

The theory has been implemented with some success in what may be 

called multiplicative systems. If x, the signal to be processed, 

cc 3 

consists of components x^, and combined as in: x = x^, * x^, ^en the 

transformation employed might be the logarithm function (real or complex) • 
What then results is, 

log x = a 'log x^ + 3' log x 2 

so that the modified components, log x^ and log Xg, have been combined 
in the ordinary linear sense. Whether or not anything has been gained, 
as well as the form of the processing to follow, depends upon certain 
properies of x^ and x^ 9 Another class of systems to which the theory 
has been applied may be called convolutional. 

If x, the signal to be processed, consists of components x^ and 
x 2 which are combined by a convolution, denoted x = x^ * x^ 9 the 
generalized linear filtering technique may be effective. The trans- 
formation might consist of, first, the evaluation of the Fourier 
Transform since it maps a convolution to a multiplication: 



F [x] - F [xj -F [xj 
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As before, this can now be followed by a logarithm function evaluation. 
Here, the logarithm must, in general, be complex. V/hat results is, 

log F [x J = log F [xj + log F [\ 2 J 

Again, whether or not anything has been gained depends on the 

4 

properties of x^ and x 2 » Schafer has examined a specific class 
of signals in which the x^ and x 2 exhibit favorable properties. His 
presentation is closely followed herein. The signals may be described 
as echoed, i.e., of the form 

k 

(1.1) x(t) = E a s(t-T ) 

i=0 1 

Such a signal consists of a basic waveform with scaled, delayed versions 
of the same waveform superimposed. This is a reasonably approximate 
description of multipath signals which can arise in many applications. 
The components x^ and x^ are identified by re-writing equation (l.l) 
as a convolution of the basic waveform and an impulse train, 

x(t) = X^t) * x 2 (t) 

(1.2) x^t) = s(t) 

> k 

x (t) = E a 6(t-T. ) 

^ i=0 

where 6(*) represents the Dirac delta function 



The important aspect of the identification made in equation (1.2) 



-26- 



is that 



log f[x^ 



is a slowly varying function of frequency while 



log F j^x 2 j is, in many cases, a rapidly varying function. Consequently, 



the inverse Fourier Transforms of log F 



L X 1 



and 



log F [ : 



are, in 



log F J^> 



operation F 



many cases, somewhat disjoint, A more important property is that 
,-l i - - ” ' - is also an impulse train. Hence, by use of the 

log F [xj j , the convolved components are transformed 
to added components which exhibit properties that may make separation 
possible. 

A similar operation was described by Bogert, Healy and Tukey ^ 
in which the power spectrum of the log of the power spectrum of the 
convolved signals was used. The result was called the cepstrum of 



the signal, denoted x, 



A 

x = 



.-1 



log 






Many properties of both operations are similar. In particular, 
the impulse train is mapped to another impulse train in both methods. 
A major difference is that, with the cepstrum, the loss of the phase 
information precludes any inverse transformation. Although filtering 
and recovery are not possible, some parameters of the signal are 
measurable in the cepstrum so that the technique finds many appli- 
cations. Due to the similarities involved, the transformation which 
retains the phase information has come to be known as the complex 
cepstrum while the other is sometimes referred to as the power 
cepstrum to maintain the distinction. 

On the subject of terminology; it becomes awkward to speak of 
"rapidly varying functions of frequency," as above, when the 
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temptation is to say "high frequency components - in the frequency 
domain*" In this regard, some phrases have been developed in the 
literature of this subject and will be used to avoid confusion. In 
this thesis, the term power cepstrum will be used to refer to the 
transformation of Bogert et al. The transformation developed by 
Schafer will be referred to as the complex cepstrum, or as simply 
the cepstrum. The cepstrum will be considered to be a function of 
quefrency. Thus, if the complex logarithm of the Fourier Transform 
of a signal has a component which is a rapidly varying function of 
frequency, it has a high quefrency component. Filtering operations 
on the cepstrum will be referred to as short pass (analogous to 
low pass but in the quefrency domain), long pass (analogous to high 
pass), and comb (the usual - but in the quefrency domain) operations. 

With the advances in technology leading to high speed digital 
computation and fast a/d conversion, and, with the development of 
the Cooley-Tukey algorithm, realization of the various mathematical 
operations called for by the technique has become increasingly 
practical. Since the implementation is via digital computer, the 
theory has been developed in terms of sampled data versions of the 
convolved signals. An in-depth presentation of the theory is 
contained in Schafer and only the relevant portions will be related 
herein. A brief overview of transform theory for sampled data signals 
is called for before the specifics are presented. 

A continuous - time signal, f(t), may be sampled every DT time 
units to produce a sequence, f(n) = f(n-DT). Such a sequence has a 
Z - Transform defined as, 



OO 
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F(z) = £ f(n) z n 

n=-*> 

The inverse operation is accomplished viai 

f(n) = ^ F(z) z 11 " 1 dz 

where G is some closed contour in the 
region of convergence of F(z) 

By this operation, f(n), a function of a discrete variable is 
transformed into a function of a continuous complex variable. The 
Fourier Transform of a discrete sequence is defined as the Z-Transform 
evaluated around the unit circle in the z - plane, 

F(e jw ) = E f (n) e" jwn 
. n=- p0 

The inverse operation is, 

f(n) - ■£- S F(e j ") d» 

-IT 

Both the Z-Transform and the Fourier Transform will be referred 

to in the development of the report since they allow easily accomplished 

evaluations of results. Neither, however, can be directly realized 

by a digital computer since they are functions of continuous variables. 

What is utilized is the so-called Discrete Fourier Transform (DFT) 

which is the Z-Transform evaluated at equally spaced points around 

the unit circle in the z - plane, 

N-l jZrkn 

= E f (n) e" N 
n=0 



F(k) 



and the associated inverse is 
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1 N-l ,i2rrkn 

f(n) = ~j~ E F(k) e N 
k=0 

All of these transforms may suffer from shortcomings referred to 
as aliasing and leakage. Additionally, the DFT is accompanied by 
the so-called picket fence effect. These, effects and means of 
counter-acting them are developed by Cooley et al and will be 
discussed in this report only when they present problems. The DFT 
is evaluated via the Cooley- Tukey algorithm, or the Fast Fourier 
Transform (FFT), 

In the case of interest, the Z-Transform of the sampled composite 
signal is, 

X(z) = X 1 (z)-X 2 (z) 



The logarithm is taken to obtain, 

X(z) = log X(z) = log X i (z) + log X 2 (z) 
The complex cepstrum, x(n), is to be, 



x(n) 




1 

2rrj 



$ c X(z) z 11-1 dz 



Here the complex logarithm presents a problem since it branches 

A 

at values of z for which ARG X(z) = j^nr. When this happens, X(z) 
has a discontinuous imaginary part and is therefore not analytic. 
Since the cepstrum will be computed via the DFT, it is required that 
ARG X(z)' be continuous for values of z along the unit circle. The 
problem is overcome via a so-called phase unwrapping technique which 
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identifies the proper branch of the arctangent function defining 
the phase. This technique is described in Appendix B. 

More specifics of the computation will be presented in subsequent 
chapters. The important aspects of the technique can be summarized: 

1, Convolutions are mapped to additions. 

2, The mapping is continuous; the inverse mapping exists 
and is continuous - hence th'e term homomorphic 
deconvolution, 

3, Deconvolution is accomplished via a subtraction (if 
it is known what to subtract) . 

4, The components of the original signal, in some cases, 
become somewhat distinctive in the complex cepstrum 
so that the may be identified and separated. 

The technique has found many applications since it can deconvolve 
signals when the components are unknown. It is merely necessary that the 
signal be reasonably echoed in nature with reasonably spaced echo 
arrival times. The basic waveform, either as it is or with minor 
processing, will have a low quefrency nature and decrease in 
magnitude at least as fast as l/n in the cepstrum. The impulse train 
can be made to have its lowest quefrency component displaced from the 
zero quefrency line by a value corresponding to the first echo arrival 
time. In cases of large echo times, the echoes can be removed by 
a simple short pass or more involved comb filter. Conversely, the 
basic waveform can be removed and the impulse train recovered by 
use, of a long pass filter. This approach is applied in cases in 
which the multipath response, x is not strictly impulsive and it 
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is desired to see precisely what it is. In such a fashion, properties 

of the paths involved can be measured. 

Some applications of the technique have been to speech proces- 

3 

sing by Schafer, Oppenheim and Stockham , to seismology, where it is 

7 

desired to recover and measure x^, by Ulrich , and to processing of 

g 

measured brain waves evoked by visual stimuli by Senmoto and Childers . 

In this last application, the data is available via an electro- 
encephalogram. There is some degree of control in that the echoed 
signal measured represents the composite response of the brain to a 
pulse-like visual stimulus. The stimulus can be re-applied, at 
fixed intervals, over a pulse repitition interval. A somewhat 
involved technique is used to average the responses obtained over 
a repitition interval so that improved signal to noise ratio is 
obtained before computing the cepstrum. 

The similarities between this application and the Loran-C 
case lead to speculation on signal to noise ratio improvement 
which will be discussed in Chapter 6. Considerations such as the 
effects of interference, atmospheric noise, and signal stability 
must be taken into account. Before these are considered however, 
more specifics of the problem must be developed. 



Chapter 2 

PROBLEM FORMULATION 
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In order to apply the homomorphic deconvolution technique, a 
model of the Loran-C signal as described in Chapter 1 must be 
developed. From this model the signal processing tasks can be identi- 
fied so that an appropriate solution scheme results. To this end, 
the waveform to be processed will be considered to be obtained 
from one pulse duration of the received signal, represented as 
f ollows : 

k 

(2.1) F(t) = n(t) + 2^ A f[ s i (t) sin [^(t-T^ ^ 

In equation (2.1), the zero subscripts identify parameters to 

be associated with the received groundwave. The subscripts 1 through 

k identify parameters of the first through k skywaves respectively. 

*th 

) denotes the normalized pulse envelope of the i component 
of the composite signal. The T^'s and A^ f s refer to the envelope 
arrival times and magnitudes, respectively. 0^ represents the 
phase shift between the envelope and the carrier of the i^ 1 
component. 

With a known carrier frequency, an equivalent method of speci- 
fying the signal is by use of the complex envelope representation 
as in equation (2.2). 

-s i (t) 



( 2 . 2 ) 



f(t) - n(t) 



£ A. exp 
i=0 1 



■> [»c T i 
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Figure 2-1 Generation of Complex Envelope 



A common method of describing how the complex envelope is 
obtained is shown in figure 2-1. An analysis of the signal flow 
in the figure shows that, 



'.w 


= n(t) cos w t - 
c 


E A, 
i=0 


s i (t) sin J 


w t. + e. 
c 1 1 


?.<*> 


= n(t) sin w t + 


k 

E A. 
i=0 1 


s^(t) cos 


r w t. + 0. 

[ c 1 1 



If the following definitions are made: 



n(t) = n(t) sin w c t - j n(t) cos w Q t 

and 



f(t) = f s (t) - j f c (t) 



it follows that, 



(2.3) f(t) = n(t) + E A, exp 

i =0 



w T. +0. 

c i i 



* ± (t) 



- 3 ^ 



which is identical to equation ( 2.2 ). Clearly, the complex envelope 
representation for continuous signals is no more than a mathematical 
convenience. In a sampled data form however, it can "be actually 
implemented and provides a compact means of handling the quadrature 
components of the baseband signal. The realization method is as 
depicted in figure 2-2. 




Figure 2-2 Realization of Sampled Complex Envelope 



The component f (t) is sampled and stored in the real part of 
s 

the sample word while f (t) is sampled, complemented, and stored in 
the imaginary part. What results is a sequence of the following 
form, 



f(n) - 



n(n) 



k 

£ 

i=0 



A. exp 



j jw 



T. 

c 1 



0 i 



(n) 



For the remainder of the presentation in this chapter, the signal 
will he considered to be noise free. The theory of noisy signals 
is developed in Chapter Consequently, the signal of interest is, 



- 35 - 



f(t) = 



k 

2 

i=0 



exp 







- 


j 


w t 4 + e. 
c i i 




- 




J 



^(t) 



For notational convenience, and without loss of generality, the 
magnitudes will be normalized with respect to Aq and the phase 
notation will be simplified: 



A' 

A i 



A i 

— . , i - 0, 1, . . . , k 

0 



*1 - “c T l + 6 i 



B. = A. 

l l 



r 



Hence, 

(2.4) f(t) 



k 

2 

i=0 



B i s i« 



In Chapter 4, distorted signals will be considered. Until then, 
it will be assumed that s^(t) = s(t-T\). In this case, equation 

(2.4) becomes 

k 

f(t) = 2 B s(t-T ) 

i=0 



This is the form of signals for which the homomorphic deconvolu- 
tion technique has been developed. Additionally, it is possible to 
identify the tasks to be performed at this point. As presented in 
Chapter 1, it is desirable to measure $ q 9 the groundwave parameter 
from which time measurements are made, and T^, to measure the difference 
in propagation speeds of the envelope and carrier. Additionally, 
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measurements of and B^ would be helpful. Finally, measurements 
of all 0^'s, 1^'s, and B^ 's would be desirable since a total 
power receiver might then be feasible. 

V/ith the problem thus defined, solution techniques can be 
considered. To begin, consider the following to be a reasonable 
approximation to the received Loran-C pulse envelope: 



s(t) = t^ e"" 0 ^ u ^(t) ; a = ) x 10^ sec ^ 



^6T J 



where u ^(t) is the Heaviside unit step function 



Initially, make the following assumptions: 



*o - 0 



T i = n i’ DT 



m = an integer 



DT = the sampling time 



The sampled waveform in this case is: 

k k 

f(n) = £ B. s(n-n. ) = s(n) <£> £ B. 6(n-n.) 

i=0 1 1 i=0 1 1 

where <g) denotes a discrete convolution and 6(n-n^) is a unit sample 



n = n. 



6(n-n^) = 






0 



n 
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Therefore, 



where 



f(n) = s(n) ® p(n) 



k 

p(n) = E B. 6 (n-n. ) 
i=0 

In the case of k = 1, i.e., in the case of a signal consisting of 
tiie ground wave and one skywave, 

p(n) = <5^)+ 6 (n-n^) 



If |b^ <1, the complex cepstrum of p(n) is (see Appendix A for the 

details of the computation) s 



OO 3^ 

$(n) = I (-l) m+1 ~ ■ i(n-mn 1 ) 

m=l 



whereas the cepstrum of s(n) is, 

(2.5) 



^(n) = -a 6 (n) + E ^ - + e -0,1 * 1 6(n-m) 

m=l 



m 



Since the complex cepstrum operation has the property of mapping 
convolutions in the time domain to additions in the quefrency domain, 
the composite cepstrum consists of s(n) + p(n). In this case, a 
simple estimator can be made to scan the cepstrum sequence and 
determine n^ from the location of the components of the impulse train. 
A simple comb filter could then be employed to remove all traces of 
the echo from the cepstrum - only slightly effecting the component 
due to s(n). Before filtering p(n) , an estimate of |B^| and 
could be made from the components of the impulse train. In such a 
manner, estimates of T^, | j 9 an< } are obtained and all skywave 
contamination is removed. 
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When is not an integer number of sampling times greater than 

T Q , f>(n) no longer has such a simple form. The components of the 

sin x 

impulse train become — - — terms which spread throughout the cepstrum 
sequence as in equation (2.6) (See Appendix A). 



( 2 . 6 ) 



p(n) 



OO 

2 (- 1 ) 

m=l 



m+l ^1 

m 



sin 


IT [n-m(n 1 +A 1 )j 


TT 


1 

3 

H k 
£ 
\ ; 



where = an integer 

and = DT*(n^+A^) 

- 0.5 < < 0.5 

In this case the filtering becomes less straightforward, A simple 
comb filter can not be effective. In the distortion free, noiseless 
case however, the cepstrum consists only of s(n), which is known, 
and p(n), which is known to within three parameters: B^, n^, and 
so that an estimator can be developed to identify p(n). When this 
identification is made, the deconvolution can be accomplished via 
a' subtraction. The form of the estimator is developed in Chapter 3 
and the effects of signal distortion and noise on its performance 
are presented in Chapters k and 5. 

The problem is further complicated when Tq is not an integer 
multiple of the sampling time 0 In this case, the s(n) of equation 
(2.5) has another term added to it, (Appendix A) 

A/ \ . / ^n+1 Ao 



(2.7) 



s’(n) = 



s(n) + (-1 y 



s(n) 



n 



; n / 0 



n = 0 
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The presence of this new term makes the estimation procedure described 
above suspect. Now, there are two unknown components in the 
computed cepstrum. Since the forms of the components are known 
however, a modified estimator can still be effective. 

A summary of the problem as formulated thus far can be made 
via figure 2-3» The received signal is processed to produce 
quadrature baseband components. These are sampled and digitally 
combined to produce the sampled complex envelope. The complex 
cepstrum is computed. From the cepstrum, signal parameters are estimated. 
One of the uses of the estimated parameters is shown: the ground wave 
phase estimate provides feedback to the local oscillator to permit 
phase tracking. An important aspect of this particular use is that 
it provides a method of comparing the performance of the homomorphic 
filter with other receiver techniques. 



ANALOG , DIGITAL 
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Figure 2-3 Homomorphic Filter Input-Output 



Chapter 3 

DECONVOLUTION OF NOISELESS SIGNALS 



-in- 



ks developed in Chapter 2, it is desired to estimate waveform 
parameters from the cepstrum of the received signal so that the 
separation, or deconvolution, can be achieved. In the distortion 
free, noiseless case, the technique is as presented in this Chapter. 

It is presumed that, with one skywave, the received signal will 
have a cepstrum of the following form: 



(3.1) 



f(n) = s'(n) + p(n) 

= s(n) + g(n) + p(n) 



where, from equations (2.5)* (2.6), and (2.7), 



m+1 



(a) s(n) = -a 6(n) + £ e - ^' 1 6(n-m) 

m=l 



m 



(3.2) 



00 s(n) 

(c) p(n) 



0 ; n = 0 

/ ^sii+1 AO ; n / 0 

' ' n 



Z (-l) m+1 


B™ 

1 sin 


'n[n-m(n 1 -»A 1 )j 


m=l 


la ■ 

TT 


_n-m(n 1 +A 1 )j 



In the estimation procedure it is desired to use the knowledge of 
the form of the above sequences in processing the received signal. 
Consequently, it will be necessary in the development to refer to the 
above sequences and corresponding components of the received signal 
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simultaneously. For clarity of presentation therefore, the notations 
s ( n )» g(n), p( n )» and f(n) will he used to refer to the predicted, or 
ideal, sequences as defined above whereas the received signal will 
be described as, 



r(n) = s (n) ® g.(n) ® p (n) 

(3.3) a a a 

or 

£(n) = s (n) + l(n) + $» 



The subscript a signifies actual components of the signal as 
received. In the distortion free case, the only difference between 
the components of f(n) in equation (3*2) and r(n) in equation (3.3) 
is due to computation error, 

A third set of sequences will also be referred to. Estimates of 
signal parameters Aq, A^» n^, and will be made from the received 
signal. From these, estimated sequences will be generated. The 
estimated sequences will have the form of those in equation (3.2), 
but with parameters as estimated from the sequences of equation (3«3)° 



A 



g(n) - g(n) 



where Aq is the estimate of Aq obtained from r(n). 






>*- 

(-D n+i 


Ao j . 

; n f 0 

n 1 


II 

o 

L< 

11 

o 

<\ 


0 


; n = 0 



3.1 Sequence of Estimation 

The received signal consists of the terms s (n), g (n) f and p (n) 

a a a 

which are combined to produce the cepstrum sequence shown in figure 

3 - 1 . 



- 43 - 




(ni+Ai) 



Figure 3-1 Cepstrum of Simulated Received Signal 
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It can be seen that there is insufficient seperation between the 
sequences to allow for precise estimation of the various parameters 
directly* Consequently* the estimation is accomplished in a number 
of steps; systematically estimating the parameters one at a time* 

The procedure used herein is to make the estimates in the following 
order. 



1. £(n) 

2. e;(n) 

3. p(n) 



an estimate of the received signal basic 
waveform. 

an estimate of the small linear phase term 
not removed before computing the cepstrum 

(a) n^ - a coarse measurement of the 
skywave arrival time 
(4) A-£ - a precise measurement of the 

skywave arrival time 

/ \ ^ 

(c) - skywave magnitude and phase 



The procedure begins by utilizing the fact that a good a priori 

estimate of s (n) is available, to wit: s(n) of equation (3.2). 
a 

Therefore the first estimate is s(n) = s(n). This estimated 
sequence can be subtracted from the received cepstrum to generate a 
new sequence denoted x(n;, 

x(n) = r(n) - s(n) = g a (n) + p a (n) + £ s (n) 

where 

6 (n) = s (n) - 5/ x 

s v ' a v ' (n; 
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In this chapter it can be assumed that this error sequence, £ (n) , 

s 

is negligibly small. The procedure as thus far discussed was accomplished 
on a simulated received signal and the resulting x(n) is shown in 
figure 3-2. 




Figure 3-2 Simulated x(n) Sequence 



It can be seen that there is improvement over the situation shown in 
figure 3-1, but that the components of g (n) and p (n) still interfere 

cl cl 

so that the desired estimates are still not readily obtained* More 
precisely, it can be seen that for quefrency values greater than 
( n i+Ai), the various x terms of f> a ( n ) interfere with one another 
to some extent. For positive quefrencies less than (n^ +A 1 ), there is 
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interference between g (n) and the first — term of p (n). 

a X 21 

What makes the estimation possible is the fact that for small 

negative quefrencies, p (n) is reduced enough so that it only slightly 

21 

interferes with g (n). The next step in the procedure therefore 

2L 

will be to estimate Aq from the negative quefrency portion of the 
cepstrum. 



3.2 Estimation 


of g a (n) 






Since 


g(n) 


- (-l) n+1 ; n < 0, 




it follows that 


A 0 " 








II 

1 


«a<- 2 > 






ii 


g a (-3) 




where Aq will be 


• 

• 

the estimate of Aq» As described above, it is 


true 


that, 


r(-l) 


" S a (-D + P a (-D + f s (-D * 


i a (-D 




r(-2) 


= i a (-2) + p a (-2) + e s (-2) ~ 


S a (-2) 



since the p (n) sequence is small for this range of quefrencies. The 
21 

estimate of Aq can therefore he obtained from, 



- 4 ?- 



(3.4) 



rtf 

A f 



£(-l) - £(_ 2 ) + r(-3) - r(-4) 
,111 
1 + ~r + ~r + t~ 






With Aq so estimated, the sequence g(n) can be generated and 
subtracted from x(n) to produce the new sequence y(n) s 





II 


$<n) - 


g(n) 


















0 


• 


n = 0 


where 


rJ 


9 




g(n) * 










(-D n+1 


Ao 

— ; 


n 0 



The procedure was applied to the x(n) in figure 3-2 to produce 
the y(n) of figure 3-3® 




A , \ 

Figure 3-3 Simulated y(n; Sequence 
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At this point in the procedure, since Tq = DT (hq+Aq) , the 
estimate of Tq can be made: 




where n^ is deduced in the phase unwrapping computation as described 
in Appendix C, 

The reason for the particular estimator used in equation (3*4) 
can be seen from figure 3-4. 




Figure 3-4 Comparison of g(n) and p(n) 
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Bo th sequences g (n) and p (n) decrease essentially as fast as l/n 
a a 

for negative values of n. Near the zero quefrency value, the g (n) 

3 . 

term clearly dominates the composite signal. As n becomes more 
negative however, the distinction between the two components is not 
as great. Consequently, only a few samples should be used to make 
the estimate and even among these few, the samples obtained from the 
more negative quefrencies should be less heavily weighted. This is 
precisely what is accomplished in equation (3.4). 

The next step in the procedure is to estimate the parameters 
of p(n), i.e., n i> an< l the sequence y(n). This is 

accomplished in three steps, beginning with the estimation of n^. 

3.3 Estimation of p (n) 

Q, 

The values of y(n) are examined for quefrencies corresponding 

to times between about 30 microseconds and 60 microseconds. The 

largest value of y(n) in this range of quefrencies will be assumed to 

correspond to a sample in the main lobe of the first X term in 

p (n). Consequently, the value of n at which this maximum occurs 
a 

will be the estimate of n^. The situation would typically be as 
depicted in figure 3"5» Since the values of y(n) will, in general, 
be complex, the plot in the figure corresponds to either the real 
or immaginary part of y(n). 
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n ^-2 n ^-1 n i n ^+1 n ^+2 

— =4 4 1 1 4 

"l*! 



Figure 3-5 Expected Form of p a (n) 

Once is determined, and can be estimated as follows. 
Consider the values of y(n) in the vicinity of n = n^ to be due 
only to the first — — — function in the expansion of p(n) . As an 
example of the validity of this assumption, consider the values of 
the first and second terms in the sum of equation (3.2b) for n = n^: 



PK) = 



sin TT - (n^*^)] 
B 1 ^[ni - (n^A^)] 



B^ 2 sin tt jn^ - 2(^4^)] 
2 tt ^ - 2(n i +/\ 1 )J 



+ • • • 



B^ sin ttA^ 



B^ 2 sin tt j^n^ - 2A^)J 
2tt - (n^-2^) 2 + 



The ratio of the magnitude of the second term to the magnitude of the 



first term will he, 
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| Ai| | B i| 2 | sin arAj 


Ai 


Kl 


sin 2 ttA^ 


2 -2A 1 | | B 1 1 jsin ttA 1 | 


- 2A t 


2 sin ttA^ 



Exponential weighting , as described in Appendix A, ensures that 




is 



less than 1. Additionally, the term 



sin 2 ttA^ 

2 sin itA^ 



is less than 1 so that the ratio is less than 



(3.5) 



N 

n l " 2 ^ll < 2 n i " 



for n^ > 1 
since | A-j^j < i 



Hence, for reasonably laxge values of n^, the approximation 
is a good one and allows the following: 



y(n r 2) 

y(n r i) 

(3.6) y(n t ) 

y( n 1 + i) 

y(n 1+ 2) 



B 1 sin rr [n t - 2 - (n^A^)] 


B^ sin tt(2+A 1 ) 


■n[ n l - 2 - (nj+A^J 


it(2+A 1 ) 


sin Trjn^ -1 - (n^t^)] 


B^ sin tt(1+A^) 


ir[ n l - 1 - (n 1 +A 1 )j 


^r(l+Ai ) 

X 


B^ sin it Jn^ - (n^A^)] 


B^ sin 


'"["l " ^ n i + Ai)] 


ST - 


B^ sin tt jn^ + 1 — (n^+A^jl 


B^ sin tt(1+A^) 


+ 1 - (nj+A^j 


tt(1+A 1 ) 


B^ sin tt f n^ + 2 - ( n i + A|)! 


B^ sin tt (2-tA^) 


TT^ + 2 - (rij+A^j 


1T (2+A 1 ) 



The solution can proceed in two steps: estimate A^, then B^ 



To 
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accomplish this, define the following variables: 



d, - y (n i- 2) 

1 iuqr 

yCn^l) 

d 2 “ 

yCn^l) 

d 3 ■ T^r 

iHnj+l) 
~ yC^) 



A i 

2 +A X 

~ A 1 
1 +A X 

A 1 

1 -A t 

~ A l 

2 - A ± 



It can be seen that four estimates of A^ can be obtained from 
the d^'s. The estimates obtained from d^ and d^ however, are made from 
values of y(n) which are closer to, or in, the main lobe of the first 
^ erm » They are therefore less sensitive to the effects of 
interference from other Sln X terms and the residuals of the "g (n) 

X cl 

estimate than the estimates obtained from and d^* It is, however, 
desirable to base the estimate on as many measurements as possible 
so that a compromise would be to make a final estimate based on a 
weighted sum of the four individual estimates* The procedure is: 

1. Compute the numbers d^, d^, d^, and d^ as in equation 

(3.7) 



Make four 


individual 


' 


2d x 


A 11 ~ 


i-dl 
- d 


A 12 = 


2 


i + d 2 



estimates as follows: 
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t 



13 



1 + d. 



ft 



14 



- d 4 

l - a., 



/° 

where A^. denotes the estimate of A^ based on the 
measurements which generated d^. 

3. Weight these estimates by the corresponding magnitudes 
of the measurements generating them and perform an 
appropriate scaling to arrive at the final estimate: 



<Y 

A, 






lk 



2 Kl 

k=i 1 K| 

When this estimate of A^ is obtained, equation (3*6) can be 
re-applied to estimate B^. To do this, make the following individual 
estimates. 



B 



11 



/A. 

Tr(2 + 1) A r 

sin\.(2 Tfj >'( n r 2 > 



"B 



12 



A 

B, 



Ik 



- '” (1+ K ) Kn,-1) 

_ ♦ „ — { A A \ J- 



13 



14 



15 



sin tt( 1 + A t ) 

fti 



sin tt. 



A 



y(nj 



rJ 



it(1 - A 

sin ir(l - A t ) y(n l +1) 

tt(2 - 2^) A 

sin it(2 - & t ) y(n l +2) 
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These are combined to produce the final estimate: 



2 

k=l 



B 1 = 5 



y(n x + 3 - k) | B 



lk 



2 

k=l 



y(n 1 + 3 - k) 






With the estimates of n^, A^» and B^, the estimated sequence 
p(n) can be generated. This and g(n) can be subtracted from the 

A ^ 

received signal cepstrum, r(n), to produce s(n). If the inverse 
DFT of s(n) is computed, what would result would be an estimate of 
the complex log of the groundwave. From the zero frequency point 
in this sequence can be measured - the phase of the groundwave. 

If there are more than one skywaves present in the received 



waveform, $ (n) will have additional 



sm x 



terms centered at 



quefrencies of » 2 ^ n 2 + ^2^’ • • • » (n^+A^) . » etc, and 

at various combinations of (n^-fcA^) and (n^+A^) as described in 

Appendix A, Once the terms due soley to the first skywave are 

sin x 

removed however, the lowest quefrency — — — term that remains is 
due only to the second skywave so that the above estimation procedure 
can be re-applied to estimate n^, an< ^ ®2* This can be done as 
often as necessary. The entire estimator can be described as in 
figure 3-6. 
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Figure 3-6 Homomorphic Filter -Estimator 
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J.h Performance of Estimators 

Before the results of applications of this technique axe pre- 
sented, a reason for the ad hoc nature of the estimators can be 
given. Due to the sequential nature of the estimation scheme and 
the manner in which the received signal is mapped to the cepstrum, 
the estimators are called upon to act in the presence of somewhat 
structured interference. As will be discussed, the interference 
structure can be anticipated so that the estimators can be designed 
to provide some degree of protection against it. To show this, it 
should not be necessary to work with the complicated estimators as 
developed. Instead, a similar, but much simpler, estimation 
situation can be considered. This will allow a clearer presentation 
of the important property of the estimators. 

As a rough model of the sequence to which the and B^ estimators 
are applied, consider the sequence z(n) shown in figure 3-7» 




Figure 3-7 Simplified Signal to be Estimated 
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The sequence is known to be of the form shown, but the value of a 

is unknown and to be estimated. The important manner in which 

z(n) is similar to p(n) is that the signs of the z(n) sequence 

values have a particular code to them: the values of z(n) axe positive 

for odd negative values of n and even positive values of n. This 

code is similar to that encountered in p(n) since, in p(n), there 

* s in X 

are two samples in the main lobe of the first — - — term, causing 
a similar sign reversal. It is desired to consider the effects of 
the presence of structured interference on the estimation of a. 

To do this, consider first, the estimator of equation (3*8) which 
is of the form used for and B^. In effect, the estimators are 
matched to, or instep with, the code of the z(n) sign reversals. 



(3.8) 



\ = -2 zJ-2) 

*2 " Sa^) 

a 3 = - . a (l) 

% = 2 z a (2) 



Now suppose that the actual sequence, z (n) , consists of 

a 

the anticipated z(n) as in figure 3“7» plus a constant interference 
sequence, w(n) = b. The sign code of w(n) , for n = - 2, -1, 1, 2, 
i& + + + +. This can be described as being orthogonal to the sign 
code of the estimators, which is + - + -• 
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In this case, z (n) 
' a v ' 



= z(n) + w(n) = z(n) + b, so that, 



& 



-2 ( - 



"I + t >> 



= a - 2b 



A> 

a 2 = 



( a + b) 



= a + b 



ro 



- ( -a + b) = a - b 



a^ = 2 ( 2 " + t) = a + 2b 



Hence the individual estimates have complementary error. It should 
be assumed that j a j >2|b| and a> 0 in making the following 
argument about the combined estimate. 



^ a (-2) 


1 A-* 

r a l + 


i 2 

1 a 


(-D | 


I ^ . 

|- a 2 + 


| Z a (1) 




\ + | ! 


>a< 2 >i 


' a 4 


1 


z a<- 2 >| 


-f 


Iz (-1 
a 




1 + 1 


z a"h[ 


+ 




r 



Since a > |2b|, it follows that, 



l (- 2 )| ■ |( - -|- + ») 



(-1) = | (a + b)| 



z (1) 
a N ' 



= |(a - b) 



z (2) 
a x ' 



(-F + ») 



(a + b) 

(a - b) * 

<-r + b > 
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From this, it also follows that 



rJ 

a 



(~ - b) (a - 2b) + (a + b) 2 + (a - b) 2 + (^ + b) (a + 2b) 



(§ - b) + (a + b) + (a - b) + (| + b) 



(3.9) 



a 



3a 2 -I- 6b 2 
3a 




As an example, suppose that b = a/5: 



r-> 

a 



a + ~~ = 1.08 a 



25 a 



The estimate for this technique is 0.08a whereas if only z (1), 

a 

for example, were used, the error would have been b = 0.2a, If b = .la, 
the errors are 0,02a and 0.1a respectively. If it happens that the 
interference has the same code as the estimator, it can be shown 
that the error is 0.38a for b = ,2a and 0.17a for b = .la. 

If w(n) = (-l) n b, the individual estimates will be, 



a^ = -2(- ~~ + b) 



a - 2b 



a, 



'2 



(a - b) 



a, 



= - (-a - b) 



a + b 




2 < -f- + »> 



a + 2b 



Once again the errors are complementary and, it can be shown, the 
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total estimate will be as in equation (3.9) • Figure 3-8 shows the 

variation of the estimate error as a function of b for an interference 

sequence of constant magnitude. The three plots are of; estimate 

from z (l) alone, four step estimate - interference in step with 
a 

estimators, four step estimate - interference orthogonal to 



estimator. 




Figure 3-8 Estimator Performance Curve 
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In the actual case, the estimators axe a bit more complicated 

than those of equation (3 .7) and the symmetry is not quite complementary 

as described. The example does however, illustrate the improvement 

obtained by centering the measurements around the main lobe of the 
sin x 

— - — functions. Two cases wherein this property is helpful can 
be presented. 

The use of the estimated sequence g(n) will generate some 
small residual: 



£ (n) = g (n) - g(n) 

o U 



^O'V (.!)»+! „ (.!) 



n 



~ a 0 

n 



This is one source of interference when the and estimates are made. 

Although the sequence is not of constant magnitude, the change is 

only fron (^ Q )/6 to (^a q )/10 in the quefrency range of interest 

for n^ = 8. This is a reasonable enough approximation to a constant 

magnitude to produce some improvement. A more important case, since 

the magnitudes involved will typically be greater, is the interference 

from the second — — term. As developed in equation (3»5)t the 

magnitude of this term will be fairly well reduced for quefrencies 

near n^. The situation can be examined a little more closely to 

s in x 

see that, if the largest value of the first — — — term is denoted a, 

where B 1 sin 

a = 9 

TT 

then the interference sequence, denoted w(n) , due to the second 
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sin x 



function has values, 



| w (n 1 +2) 


1 = 


| w(n^+l) 


= 


• 




• 




1 • 

|w(n 1 -2) 


as 


CO 

• 

0 

1! 


n l = 



B^ 2 sin ttA^ 
2v(n^ + - 2) 

B^ 2 sin ttA^ 
ZnCn^’ + Ajl - 1) 



B^ sinnA^ 
2rr(n^ + A^ + 2) 



= a B. 



= a 



B, 



'll 2(n J _ + A 1 - 2) 

2(n 1 + A 1 - 1) 



= a I B. 



2 + A x + 2) 



0.031a to about 0.019a. Once again, the symmetry is not perfect, 
but it is close enough so that some improvement can be gained. 



3.5 Application to Simulated Signals 

In the case of the generated received signal used for the 
example depicted in figures 3-1 through 3-3 » the following parameters 
were used: 

DT = 5 (microseconds) 

T 0 = ll«5 

T 1 = T 0 + 39,5 J n l = 8, A 1 - -0.1 

B 1 = 0.8 

^0 = 0 



When the technique developed herein was applied to this signal, the 



resulting estimates were: 



11.48785 

8 
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X 1 

K 



-0.096033 



'= T q + 39.5198 
= (0.803417, -0.000881) 
= 0.0006864 radians 



0.803417 - jO. 000881 



The phase error in this case, \$q ~ ~ $Qt 1® 0.039 in 

degrees. At a carrier frequency of 100 kHz, this corresponds to a 
time error of about 1.1 nanoseconds. The procedure was repeated 
for a variety of cases and the results are presented in Table 3~1» 



*0 


T 

0 


T 1 


B 1 


rj 

T o 


/N-> 

T 1 


B 1 


(rad) 


error 

nanoseconds 


0 


10.0 


40.0 


.8 


10.0004 


39.9996 


*‘-3- 
O O 
O O 
O O 
O O 
O O 
CO - 
• 


-.000043 


-.068 


0 


10.0 


39.5 


.8 


9.9404 


39.4404 


(.810026, 

-.00015) 


.000077 


.044 


0 


10.0 


38o0 


.8 


9. 8091 


37.9863 


(.811575, 

-.00036) 


.000251 


.400 


0 


11.5 


40.0 


.8 


11.5326 


39.9815 


(.79 13W, 
-.00017) 


.000125 


.195 


0 


11.5 


39.5 


.8 


11.4878 


39.5198 


(.803417, 

-.00088) 


.000686 


1.095 


0 


11.5 


38.0 


.8 


11.3366 


37.9698 


(.803411, 

,00019) 


-.000227 


-.361 


0 


8.5 


40.0 


.8 


8,5685 


39.9902 


(.795314, 

.00024) 


-.000239 


, 

O 

CO 

• 


0 


8.5 


39.5 


.8 


8.5189 


39.5029 


(.799695, 

.00054) 


-.000464 


-.740 


0 


8.5 


38.0 


.8 


8.3650 


37.9386 


(.808735, 

.00013) 


-.001570 


-2.50 



Table 3-1 Results of Simulation: No Distortion, One Echo, Zero Phase 



(Note: For convenience, the T^ column in the tables actually represents 
skywave arrival time relative to the groundwave; i.e,, (Table T^) = T^-Tq) 
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The results shovm in the table reflect the performance of the 
estimators on signals which were purely real to begin with. The 
phase errors would therefore define the precision available in the 
loop shown in figure 2-3 in which the homomorphic filter-estimator 
serves as the phase sensitive detector. 

The results of the simulation for non-zero phase cases are 
shovm in Table 3“2. What can be seen in these cases is that the 
B^'s are essentially -.19723 radians out of phase with the true values,. 
This can be explained by examining the Fourier Transform of p(n) : 

P(e^ w ) = B Q + e“^ vm l = B Q (1 + e~^ m l) 



In other words, all magnitudes and phases in the complex cepstrum 
are measured relative to the groundwave magnitude and phase. 



^0 


T 0 


T 1 


B 1 


A J 

T 

A o 


T 1 


B i 


? 0 ( rad ) 


error 

nanosec 


.19723 


8.5 


40.0 


.8 


8.5687 


39.9900 


(.779802, 

-.155877) 


. 196285 


-1.50 


.19723 


8.5 


39.5 


.8 


8.5169 


39.5030 


(.783872, 

-.157803) 


.197302 


.112 


.19723 


8.5 


38.0 


.8 


8.3671 


37.9392 


(.792381, 

-.160622) 


.202399 


8.08 



Table 3-2 Results of Simulation: No Distortion, One Echo, Non-zero Phase 



The cases in which the skywaves are larger than the groundwave 
are considered next. As discussed in Appendix B, exponential weighting 
must be used in these cases. The only thing that changes in the 

AJ 

estimation procedure is that &(n), the a priori estimate, must also 



show the effects of the exponential weighting. This is easily 
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accomplished as shown in Appendix A. The results are presented in 
Table 3-3. The extra column, BJ, gives the effective skywave coefficient 
relative to the groundwave after exponential weighting. 



*0 


T 0 


T 1 


B 1 


B i 


T o 


T 1 


A J 

B i 


$ 0 (rad) 


O 

CD 

3 


0 


10,0 


40.0 


(2.0, 

2.0) 


(.51332, 

.51332) 


10.0001 


39.9999 


(.513279, 

.513268) 


.000024 


.039 


0 


10,0 


39.5 


(2.0, 

2.0) 


(.52213, 

.52213) 


9.9419 


39.4499 


(.523293, 

.533561) 


-.00581 


-9.30 


0 


10,0 


38.0 


(2.0, 

2.0) 


(.54945, 

.54945) 


9.8493 


38.0031 


(.547588, 

.560679) 


-.00357 


-5.7 


.78540 


10.0 


40.0 


(2.0, 

2.0) 


.725 


10.0000 


39.9841 


(.719902, 

.001483) 


. 784204 


-1.90 


.78540 


10.0 


39.5 


(2.0, 

, 2.0) 


.74 


9.9450 


39.4468 


(.748461, 

.001661) 


.784109 


-2.0 5 


.78540 


10.0 


38.0 


(2.0, 

2.0) 


.775 


9.8156 


37.9970 


(.785317, 

.000722) 


. 784751 


-1.02 



Table 3-3 Results of Simulation: No Distortion, One Echo, Skywave 
Larger than Groundwave 



In the cases presented thus far, either the groundwave and the 
skywave were in phase, or one of the two was at zero phase. For 
completeness, Table 3-4 presents the results for different, non-zero 
phases o 



. 


T 

0 


T 1 


B 1 




kJ 

T 0 


A J 

T 1 


\ 


$0 (rad) 


I0q error 


.19723 


10.0 


40.0 


(£.-£) 


(.39233, 

-.58828) 


10.0002 


39.9993 


(.392375, 
--58821 1 


.197271 


.064 


.19723 


10.0 


39.5 


Tf^¥) 


(.39233, 

-.58828) 


9.9412 


39.4472 


TT388902, 

-.59842) 


.201864 


7.36 


,19723 


10.0 


38.0 


(1,-i) 


(.39233 

-.58828) 


9. 8134 


37.9994 


(.38774, 

-.59804) 


. 199813 


4.10 



Table 3-4 Results of Simulation: No Distortion, One Echo, Skywave 
and Groundwave at Arbitrary Phases. 
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In all of the above cases, the final estimate of the groundwave 

A J 

cepstrum sequence, s(n), was short pass filtered after 150 micro- 
seconds before taking the inverse DFT to measure 0q» As discussed 
in Appendix A, this is done to facilitate filtering in the multiple 
sky wave case. Table 3~5 shows the results of the estimation scheme 
when applied to a signal with two skywaves. The results axe not 
quite as good as in the single skywave ca'se. Nevertheless, it is 
felt that, if warranted, more elegant estimation algorithms can be 
developed. The point of the presentation thus fax is to demonstrate 
that, in the distortion free, noiseless case, measurement precision 
on the order of a few nanoseconds is possible. The effects of 
signal distortion can now be examined. 
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Chapter 4 
DISTORTED SIGNAIS 
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The results of the proceeding chapter depend heavily on the 
assumption that the received signal is identical in form to the 
sequence f(n) defined in equation (3.1)s 

f(n) = s(n) + g(n) + p(n) 

In practice, it can be expected that g (n) = g(n) since this component 

cl 

is merely due to the small portion of the linear phase not removed 

before computing the cepstrum. There remain however, two sources 

of error since s (n) and p (n) may not be of the expected form. The 
a a 

sensitivity of the estimators to these errors will be examined in 
this chapter, 

4,1 Distorted Basic Waveform 



In the estimation scheme of Chapter 3 it was assumed that a good 

a priori estimate of s (n) was available. What is meant by good might 

a 

be that some criterion such as equation (4.1) is met. 



(4.1) 



c.W 




s a (n) - s(n) 


FTiT 







for all n such 

that s (n) ^ 0 
a 



where £ is some appropriate bound. Upon close examination however, it 

can be seen that what is of importance is the effect of the residual 

sequence, £ (n) , on the sequence x(n) • 
s 
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- & a (n) + g a (n) + A p a (n) - I(n) 



= g a (n) + p a (n) + 6 s (n) 



Consequently, what is of concern is the ratio. 



(4.2) 



e x (n) e S M 

l a (n) + jp a (n) " 6 a ( n ) + P a ( n ) 



for all n of 
interest 



As an example of the distinction, the ratio of equation (4.1) 
might have a value of about 0.01 for some value of n. As this stands 

A J 

by itself, £(n) might therefore be considered a reasonably good 
estimate. However, if s (n) is 100 times greater than p (n) + g (n) 

3, cl a 

at this value of n, the ratio of equation (4.2) would equal 1. The 
estimate might then not be considered so good - particularly if the 
value of n in question happened to be near n^. 

Conversely, if the ratio of -equation (4.1) is quite large for 
some quefrency value not considered in the estimation procedure, or, 
for some quefrency at which p (n) + g (n) is much larger than s (n) , 

ci EL 3i 

there may be no problem,, The criteria for judging the appropriateness 
of the a priori estimate would therefore be as in equation (^#3)» 

n = -1, -2, -3, -4 

n x n^ 



( 4 . 3 ) 



£ s ( n ) « g a (n) 
£ s (n) « $ a (n) 



Distortion of s (n) will be 
a ' 



introduced by the propagation path 



and by the receiver processing. It may be assumed that the propagation 
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distortion of the ground wave is negligible. Not so negligible 
however, would be the effects of receiver filters. In addition to 
performing the low pass filtering in the demodulation scheme shown 
in figure 2-2, the filters are called upon to suppress atmospheric 
noise and interference from other navigation or communication systems 
operating in the low frequency portion of the spectrum. Since no 
general linear time invariant filter can be used to model the effects 
of these filters under all conditions of operation, what will be 
presented will be a method of analyzing the technique which could 
then be applied to any specific implementation. 

Suppose the received signal is processed by a single pole low 
pass filter (similar arguments can be applied for double or triple 
pole filters). A model of the basic waveform, i.e., the signal 
without echoes, to be sampled would then be, 



s a (t) = s(t) * h(t) 



where 



h(t) = e ^ u ^ (t) 



so that 



\(n) = s(n) + h(n) 



and 




In this case, the cepstrum of the filter response is easy to 




n=0 



[l - e" ( ^ w) 



compute : 



-1 



oo 
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and 



log H(e JW ) 



2 

n=l 



e 



-Pn 



n 



e -jun 



for p > 0 



£ s (n) - h(n) 



0 ; otherwise 




; n > 0 



The fact that the cepstrum is zero for negative quefrencies will 
be true for any minimum phase filter - any linear phase term being 
removed in the phase unwrapping procedure. Since the Aq estimate is 
made from the negative quefrency values of the received signal, it 
can be seen that the accuracy of the estimated sequence g (n) will 

cl 

not be affected by the presence of the h(n) term. The values of 

A 

C (n) for quefrencies near n 1 however, will cause problems. The 

S X 

values of 6 (n) and p(n) in this region are given in Table 4-1 for 
s 

comparison. The parameters of the signals in this case are: 



P = a - 2/65 

n^ =8 

1 = - 0.1 

B 1 = 0.8 

DT = 5 



n 


a — ; — r 

e K ( n ) 


p(n) 


n.-2 


.13857 


-.0348 


-n r l 


.11518 


.0792 


n i 


.09770 


.7796 


n^+1 


.08423 


-.0791 


n.j+2 


.07351 


.0472 



Table 4-1 Comparison of A Priori 

Estimate Error and p(n) . 
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Cleaxly, 



€ s (n) 



is too large in this case to permit estimation 
of A^ a nd. B^. The problem can be overcome if the response of the 
filter can be presumed to be known reasonably well. The a priori 



estimate of s (n) can then be modified: 
a 



so that 



s(n) = s(n) + h(n) 

? s (n) = s a (n) - S(n) 

= s(n) + h(n) - s(n) - h(n) 

- h(n) - h(n) 



In this case, an estimate of the form 



a / \ 
h(n) 



0 ; otherwise 

e -Yn 

; n > 0 



where y « p might be used. The error therefore would be, 



€ (n) 
s' ' 




; n > 0 



3 

For y = ~ {3 , with the other parameters as in the preceeding example, 

the error values would then be as shown in Table 4-2. 
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n 


£ s (n) 


P(n) 


n^-2 


.01222 


-.0348 


n r i 


.01176 


.0792 


n l 


.01129 


.7796 


n^+1 


.01089 


-.0791 


n^+2 


.01048 


,04?2 



Table 4-2 Comparison of Revised 
A Priori Estimate 
Error and p(n) 

A / i ^ / * 

Within this range of quefrencies, £ (n) < p (n) , Further 

S cl 

A 

notice should be taken of the fact that the € (n) terms are of 

s 

relatively constant value and of constant sign. As developed in 
Chapter 3» the individual estimates of and B^ will therefore 
have nearly complementary errors. Consequently, the total estimates 
of and are relatively unaffected by the incorrect a priori 
estimate of (3. 

If the correct value of Y is used, i.e., Y = 8, but an incorrect 
filter form is assumed, the results axe quite different. Suppose 



h(t) = t e ^ u ^(t) 



It can be shown that the resulting cepstrum is, 



h(n) 



- 




■> 


0 


; n 


< 0 


- 8 


; n 


= 0 


-8n 






2 


; n 


> 0 


n 




J 
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Therefore, 

e s ( n ) = 



This will produce the same situation as that depicted in Table 4-1, 

In summary, it can be said that the a priori estimate of s (n) 

ci 

must be changed to reflect the effects of the responses of any filters 
used. The estimation procedure of Chapter 3 will be relatively 
insensitive to errors in filter parameters, e.g, , pole locations, 
but fairly sensitive to errors in estimated filter form. 

4.2 Distorted Echoes 

If the received signal is processed by a linear time invariant 
filter, the effects will be as described above. Such a filter will 
not change the echoed property of the signal, i.e., p (n) will not 
be changed. Distortion of p (n) can only be introduced by differences 

cl 

in the propagation paths 0 The distortion can be modelled as follows 
in a single skywave case. 



0 

-e 

-£n 



; n < 0 

; n = 0 

: n > 0 



J 



r(t) = s a (t) + B i S a ( t-T l) * 



R(e JW ) - S a (e^)-G a (e JW ) 



1 + 



3 1 e-j“( n l +i L ) H(e j ") 
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log R(e J ) = log S (e JW ) + log G (e JW ) 



» v .. B k e- J " k(n l +A l ) . . 

+ £ (-J) k+1 -i j H k (e^") 

k=l 



for I H(e^ W ) | < 1 



r(n) - 



s a (n) + g» + Un) 



where 



P_( n ) 



= F 



-1 



“ B, k t . 

2 (-l) k 1 — r ll k (e JU ) 

k=l K 



The resulting p (n) can be anticipated by recognizing that its 

cl 

components are inverse Fourier Transforms of products. Hence, what 
sin x 



was the first 
h(n) = h(t) 



term in the expansion is now convolved with 



The second 



sm x 



terra is convolved with h(n) 



t=n 'DT 

twice and the third term is convolved with h(n) three times, etc. 

The results of these convolutions are combined to form p (n) 0 In 

a 

the cases considered thus far, h(t) has been an impulse function 
so that the — terms appeared intact. If h(t) is only approxi- 
mately impulsive, the first S - ^ X term will be only slightly 

sin x 

distorted. Higher order — — — terms however, will be convolved with 
h(n) a number of times and will tend to spread out - only loosely 
resembling * functions. 

In the Loran-C case, h(t) is nearly impulsive. Additionally, 



since the cepstrum is short pass filtered, the effects of the many- 

sin x 



fold convolutions which distort the high order 



terms will 




- 
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only slightly be felt. The selection of an appropriate model for 
h(t) is confounded by the fact that it is so nearly impulsive that 
it has not been possible to measure its form. Perhaps the only 
correct model would be that the skywave and groundwave are identical 
to within 1 % along the leading edges of the pulses. 

Although an appropriate form of h(t) is not available, it is 
desired to test the sensitivity of the estimators to a possible 
1 % distortion. What will be used in this case is a groundwave of 

p ry x p _Q x 

the form t e and a skywave of the form t e -9 . If a = 2/65 
and 3 = 2/70, and if the signals are scaled to achieve a value of 
1 at their respective peaks ( 2/ct and 2/3), the skywave will have a 
magnitude of about 7 of the magnitude of the groundwave at t = 65 . 
(with both signals beginning in coincidence). While the h(t) that 



r -l 



[(s-Kt)/(s-t3 )] 



3 



, is not 



would produce such distortion, i.e«, L 
claimed to be a particularly realistic model of what actually occurs, 
it does serve as a clearly large overbound on the amount of distortion 
to be expected. 

The case described above was simulated and the results are 
presented in Table 4-3 o While the errors are, understandably, 
higher than in the cases presented in Chapter 3* they are the results 
of distortion which is much more than can be expected to be introduced 
by the propagation path. As long as the propagation path filter is 
nearly impulsive, the effect on the cepstrum estimators is minimal. 

The results presented in this chapter will be discussed further 



in Chapter 6« It can be said however, that, with the assumptions made 
in the presentation, signal distortion only slightly affects the 
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performance of the homomorphic filter as used in this application. 



^0 


T 

0 


T 1 


B 1 


B i 


T 


T 1 


s i 


y > 0 rad 


>0 

error 


.19723 


10.0 


40.0 


( 2 ., 

- 2 .) 


(.32977, 

-.49448) 


10.0016 


39.875 


(.31743, 

-.47539) 


.16990 


43.5 


.19723 


10.0 


39.5 


( 2 . , 
- 2 .) 


(.33627, 

-.50421) 


9.9506 


39.4662 


(.34000, 

-.52213) 


.18663 


16.9 


.19723 


10.0 


38.0 


( 2 ., 

- 2 .) 


(.35652, 

-.53459) 


9.8317 


38.0789 


(.35804, 

-.55027) 


.18418 


20.8 



Table 4-3 Results of Simulation: Distorted Echo 



Chapter 5 
NOISY SIGNALS 
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In this chapter the effects of noise on the cepstrum estimators 
will be examined. Before quantitatively presenting the effects of 
the noise, an attempt should be made to develop a theory which 
outlines the results to be anticipated. The development begins 
with an explanation of a matrix notation which seems best suited 
to a clear presentation of the elementary aspects of the problem. 

If a continuous-time function, x(t), is sampled every DT time 
units, the resulting samples have values x(n*DT). These samples 
produce a sequence of values denoted x(n) as used in the preceeding 
chapters. The values can also be denoted x r and used to define a 
vector xs 



X = (x 0 ,x 1 ,x 2 ,...x n _ 1 ) T = (x(0),x(l),x(2 ), ...x(N-1)) T 



The sequence has a Discrete Fourier Transform, denoted x^\ 



N-l _ j2Ttkn 
x T = S x e N 

/Jk a n 

n=0 



The values of the DFT can also be represented as defining a vector 



x j 






In matrix notatation, the DFT relation can be expressed as: 
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or. 



x = W x 



( 5 . 1 ) 



r " " 
o 

x< 




1 


1 


1 


•> 

. . i 




X 1 

A/1 




1 


W 


w 2 


. . w^ 1 




x_ 




1 


w 2 


w 4 






a/2 

• 




• 


• 


• • 


• • • 


• 


• 

• 

*N-1 




• 

• 

1 


• 

>-1) 


• 

W2(N-1) 


• 

>-l) 2 

J 















"N-l 



where W = exp(- • 

It is desired to trace the effects of the noise from the time 
domain to the quefrency domain. The first step therefore is to 
examine the behavior of the DFT of a noise sequence. To facilitate 
the derivation, the noise sequence to be encountered will be considered 
to consist of samples of a zero mean, complex, white gaussian random 
process. The appropriateness of this model is discussed in Appendix D. 
To proceed, make the following definition, 



x = v ; v is a zero mean, complex, white 

gaussian random vector 

9 

The sequence x then has the following properties : 



E 



^xj = E j^Re(x) j = E (^lm(x) ] = 0 
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^ ; I_ = N x N identity matrix 

I 



where * denotes complex conjugate 

In short, x is a vector of N complex random variables. The vector 
contains a total of 2N uncorrelated zero mean gaussian random 
variables (GRV's), The uncorrelated property, for zero mean GRV's 
implies independence. The complex variables, x^, can also be represented 
as, 

X i = | x i |‘ exp (aRGCx^)] 

where |x^j is a Rayleigh distributed random variable and ARG(x^) is 
a random variable uniformly distributed on (~ o 

One way to arrive at the statistics of x is to argue along 
the following lines. From equation (5*1), 

. ir k ^ „2k ^ , ,Jc(N-l) 

x k = x 0 + V X 1 + ¥ x 2 + ... + vr x H _ 1 

(5.2) = x^ + x^ + x^ + ... + x^_ 1 

where x^ = The purpose of defining this new sequence x^ is 



E |Re(x^) Re(x) = 



_G 

2 * 



^Im(x 



E Im( xl ) Im(x) 



a 

2 



E 



ReCx 1 ) Im(x) 



= 0 



E 



T J ^2 T 

XX* = 0 I 



as follows: 
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|x^|is therefore Rayleigh distributed with statistics identical to 
those of jx^|. Furthermore, ARG(x!) = ARG(w^) +ARG(x^), Since 
the phases are computed as principal values, ARG(x^) is also uni- 
formly distributed on (-Tr,Tr). Consequently, x? has Rayleigh magnitude 
and uniform phase values which are statistically identical to those 
of x^. It follows that the real and imaginary parts of xl^ must also 
be statistically identical to those of x^, i.e., independent zero 
mean, identically distributed GRV's, It is also true that the 
elements of the x' vector are independents 



fx' x'*| = E 


„mk 
W x 


,.-nk „ 

W x* 


l m n J 


m 


nj 



= E [x x*| = a 2 6 

L m n J mn 



By the zero mean GRV property, this implies independence. The result 
of this intermediate step can be seen by re-writing equation (5. 2) 
as below. 



(5.2) 



£k = x 0 + X i +x 2 + — +X ’ N -1 



Re(x k ) = Re(x^) + Re(x’) + ... + Re(x^_ 1 ) 

(5.3) 

!m(Xk) = Im(x^) + Im(xp + ... + Im(x^_ 1 ) 
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The right sides of equation (5«3) contain a total of 2N independent 

zero mean GRV's. The two sums, i.e., R(x k ) and , are therefore 

independent. Additionally, they each are sums of K independent zero 

2 ' 

mean GRV's with variances — . It follows then that Re(x^) and 

Im(x, ) are independent zero mean, identically distributed, GRV's with 
2 

variances N — so that 

It is also possible to show that the individual elements of 
the x-vector are independent. This follows from the orthogonality 
of the rows of the W-matrix, i.e., _x^ and x^ become orthogonal 
random variables. Since they are also zero mean and gaussian, they 
are independent. To see this more explicitly, consider: 



jx^| is Rayleigh and ARGCjs^) is uniform. 



E 



(jk S] - E (*0 + + W 31 x 2 + ... + W k(B ' 1) Vl 



„ , ,.-m „ , , ..-m(N-l) „ 

x* + W x* + ... + W v ' x*_^ 



Since E fx. x*| = o 2 6. ., it follows that 

L i J J ij 



E 



x, x* 



~'m 



= E 



c Q | 2 + |x 1 | 2 + ... + w ( n - 1)(^) )Xn _| 



o N-l /T N 

0 2 E w n(k-m) 
n=0 



For k = m, the value of the above sum is clearly Ncr. For k ^ m, the 
sum consists of N terms, all with equal magnitudes and with phases 
equally spaced from -tt to tt. The resultant sum is therefore zero so 
that, 



n(k-m) 



= No 2 6 



- 8 > 



E 




x* 



~m 



N-l 

E 

n=0 



W 



km 



Once again, together with the zero mean gaussian property, this implies 
independence. 

In summary, a zero mean, complex, white gaussian noise sequence 
has a DFT which is also a zero mean, complex, vihite gaussian noise 
sequence. All that changes, statistically, is that the variances 
are changed from O 2 to Nc? 2 . 

In computing a DFT, it is sometimes desired to add zeroes to the 
end of the sequence. Since the same DT is used, this does nothing to 
change the bandwidth of the DFT. However, since there are more samples 
of the time sequence, the DFT will have more samples in it. This has 
the effect of increasing the resolution of the DFT. The effect of this 
operation on the statistics of the noise DFT should be examined. 

Suppose that M samples of the noise sequence, x n , are available 
and that an N point DFT is computed: 
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w 2 


W N-1 
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*2 
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w 2 


w 4 
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X 2 
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The x-vector will, in general, have N non-zero elements which, 

/sj 

as above, can be shown to have real and imaginary parts which are 

2 

independent, zero mean GRV's with variances M . Hence the 
magnitudes of each point are Rayleigh distributed and the phases 
are uniformly distributed. The sequence, however, is no longer 
white. In fact, the x-vector will span the same M dimensional 

/V 

complex (or 2M dimensional real) space as that generated by the 
original sequence x. The autocorrelation function of the DFT can 
be computed in a straightforward manner. From this it can be 
shown that, for example, if N/M = 4, points in the DFT 4 samples 

apart are independent. If N/M is not an integer, all the points 
axe correlated, although only loosely for seperations of more than 
about 2N/M. 

If x = £ + v, where s is a deterministic vector and v is the 
noise vector as previously described, it follows from the linearity 
of the DFT operation that 

E [2k] ‘ E [s fc + X k ] - 8* 

and 

E (x, x*| = E f(s. + v ) (s* + v*) 
l.Hk ~mj [ v '“k ~k y v ~m ~m^ 



= s. s* + No 2 6 t 
~ k km 



With this background, consider the DFT of the sequence 



* n = vfhich may appear as in figure 5-1 
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The solid line in the figure represents the magnitude of the 
DFT of the deterministic sequence. The dashed line represents some 
convenient probabalistic bound on the magnitude of the DFT of the 
noise as determined by its Rayleigh distribution. What must be con- 
sidered is the effect of the noise on the phase of the composite 
signal. For values of k such that |s^| > jvjJ , the vector sum may 
be as shown in figure 5~2(a). 




- 86 - 



In this case, the composite phase will follow that of s^. reason- 
ably closely. As an example, if |SjJ = 2|Vfc|» the phase error, the 
difference between arg(s^) and arg(s^ + v^) is at most about 31° • The 
phase curve is therefore noisy; but the noise is small - less than 
tt/2 in magnitude. 

For values of k such that |s k | < |v^|, the vector sum may be as 
shown in figure 5~2(b). The phase error in this case can be considered 
to produce a more noisy curve - errors as great as tt in magnitude. 
However, the implications of this situation are more significant in 
view of the fact that the phase curve will now have the structure of 
the noise phase, i.e. white. This causes problems because it is 
presumed in the phase unwrapping algorithm that the phase has a smooth, 
or continuous, structure. If js^j > |y k jfor all k, every point in the 
DFT phase will be within + — of the correct value so that the measured 
number of branchings of the arctangent function defining the composite 
phase will be, at most, one different from the correct value. If 
|s k j < [y k ^ however, there can be any number of branchings of the phase 
so that a huge error can be expected in the corrected phase curve. 

This can be summarized by considering the form of the complex 
log of the DFT. 



f 1 




Xk" 


lo ® [2fc + \ 


= log s k + log 


1 + — £ 



l0S Sk 








- 87 - 



los frk + ^k] 


= log X k + lo s 


- 

S T 

1 + -^L 

V. 






"k. 



log X k + 






v, 

- v 'k 



v » s, 
r k| pkl 



If the phase of v, is uniform and white, the phase of will also 

K % 



be uniform and white, when measured via its principal value. Therefore, 

Consequently, 



v T 

is statistically equivalent to 

®k 



R e (Y k ) + 3 Ira(Y k ) 



Pkl 



the effect of the log DFT operation for |s k |» |v^| is to filter the 
noise as in figure 5~3(a). The | used is that of the Loran-C signal, 




Figure 5-3 LOG DFT Filtering Effect 
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For j | < |y k |, the composite log DFT is essentially that of y k « 

Consequently, the resultant log DFT can be approximated as: 





log v, 



Me 



jARG(v k ) 



■ H, 



where and are as shown in figure 5-3 (b) and (c) and |y k | and 

ARG(y k ) are the Rayleigh magnitude and uniform phase as previously 

derived. The equivalent operation is depicted in figure 

The implications of this problem can now be seen. The output of 

the phase unwrapping block of figure 5-^ will be grossly in error if 

the Rq of figure 5-3 corresponds to a frequency less than the sampling 

frequency, 1/2BT. For a given signal to noise level therefore, sampling 

faster than a prescribed rate has dire consequences. It is, however, 

desired to sample the input signal as fast as possible, obtaining as 

many samples of the groundwave as possible before the first skywave 

arrives. The number of samples so obtained will define n^ in the 

complex cepstrum, the quefrency seperation between the terms 

of p (n) . If this number is reasonable, e.g., 4 or more, the estimators 
a 

can be expected to perform well. If it is less than 4, there will 

be too much overlap between the terms of p (n) to obtain the desired 

a 

precision'. If a greater bandwidth is necessary, some sort of 
smoothing would have to be applied to the log magnitude and phase 
curves to reduce the errors introduced by the false branchings due 
to noise. In such an operation, however, precision is lost in the 
cepstrum. 
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Figure 5-^ Development of Noisy Cepstrum 



